home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / ts_coef.pro < prev    next >
Text File  |  1997-07-08  |  3KB  |  113 lines
  1. ; $Id: ts_coef.pro,v 1.7 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1995-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       TS_COEF
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the coefficients used in a Pth order
  11. ;       autoregressive time-series forecasting/backcasting model. The 
  12. ;       result is a P-element vector whose type is identical to X.
  13. ;
  14. ; CATEGORY:
  15. ;       Statistics.
  16. ;
  17. ; CALLING SEQUENCE:
  18. ;       Result = TS_COEF(X, P)
  19. ;
  20. ; INPUTS:
  21. ;       X:    An n-element vector of type float or double containing time-
  22. ;             series samples.
  23. ;
  24. ;       P:    A scalar of type integer or long integer that specifies the
  25. ;             number of coefficients to be computed.
  26. ;
  27. ; KEYWORD PARAMETERS:
  28. ;  DOUBLE:    If set to a non-zero value, computations are done in double 
  29. ;             precision arithmetic.
  30. ;
  31. ;     MSE:    Use this keyword to specify a named variable which returns the 
  32. ;             mean square error of the Pth order autoregressive model.
  33. ;
  34. ; EXAMPLE:
  35. ;       Define an n-element vector of time-series samples.
  36. ;         x = [6.63, 6.59, 6.46, 6.49, 6.45, 6.41, 6.38, 6.26, 6.09, 5.99, $
  37. ;              5.92, 5.93, 5.83, 5.82, 5.95, 5.91, 5.81, 5.64, 5.51, 5.31, $
  38. ;              5.36, 5.17, 5.07, 4.97, 5.00, 5.01, 4.85, 4.79, 4.73, 4.76]
  39. ;
  40. ;       Compute the coefficients of a 5th order autoregressive model.
  41. ;         result = TS_COEF(x, 5)
  42. ;
  43. ;       The result should be:
  44. ;         [1.30168, -0.111783, -0.224527. 0.267629, -0.233363]
  45. ;
  46. ; REFERENCE:
  47. ;       The Analysis of Time Series, An Introduction (Fourth Edition)
  48. ;       Chapman and Hall
  49. ;       ISBN 0-412-31820-2
  50. ;
  51. ; MODIFICATION HISTORY:
  52. ;       Written by:  GGS, RSI, November 1994
  53. ;       Modified:    GGS, RSI, January 1996
  54. ;                    Added DOUBLE keyword.
  55. ;       Modified:    GGS, RSI, June 1996      
  56. ;                    Replaced nested FOR loop with vector ops. 
  57. ;                    Faster execution for values of P > 100.
  58. ;-
  59.  
  60. FUNCTION TS_Coef, x, p, MSE = MSE, Double = Double
  61.  
  62.   ;Compute the coefficients of the Pth order autoregressive model
  63.   ;used in time-series forecasting/backcasting.
  64.   ;ARcoef = ARcoef[0, 1, ... , p-1]
  65.  
  66.   ON_ERROR, 2
  67.  
  68.   TypeX = SIZE(x)
  69.   nX = TypeX[TypeX[0]+2]
  70.  
  71.   if p lt 2 or p gt nX-1 then $
  72.     MESSAGE, "p must be a scalar in the interval: [2, N_ELEMENTS(x)-1]."
  73.  
  74.   if KEYWORD_SET(Double) eq 0 then Double = 0
  75.  
  76.   MSE = TOTAL(x^2, Double = Double) / nX
  77.  
  78.   ;Do all intermediate calculations in double-precision.
  79.   ARcoef = DBLARR(p)
  80.   str1 = [0.0d, x[0:nX-2], 0.0d]
  81.   str2 = [0.0d, x[1:nX-1], 0.0d]
  82.   str3 = DBLARR(nX+1)
  83.  
  84.   for k = 1, p do begin
  85.     ARcoef[k-1] = 2.0d * TOTAL(str1[1:nX-k] * str2[1:nX-k], /DOUBLE) / $
  86.                      TOTAL(str1[1:nX-k]^2 + str2[1:nX-k]^2, /DOUBLE)
  87.     MSE = MSE * (1.0d - ARcoef[k-1]^2)
  88.  
  89.     if k gt 1 then begin
  90.       i = LINDGEN(k-1) + 1
  91.       ARcoef[i-1] = str3[i] - (ARcoef[k-1] * str3[k-i])
  92.     endif
  93.  
  94.     ;if k = p then skip the remaining computations.
  95.     if k eq p then goto, return_ARcoef
  96.  
  97.     str3[1:k] = ARcoef[0:k-1]
  98.     for j = 1, nX-k-1 do begin
  99.       str1[j] = str1[j] - str3[k] * str2[j]
  100.       str2[j] = str2[j+1] - str3[k] * str1[j+1]
  101.     endfor
  102.   endfor
  103.  
  104.   return_ARcoef:
  105.  
  106.   if TypeX[TypeX[0]+1] eq 5 and KEYWORD_SET(Double) eq 0 then $
  107.     RETURN, FLOAT(ARcoef) else $
  108.   if TypeX[TypeX[0]+1] eq 5 or KEYWORD_SET(Double) ne 0 then $
  109.     RETURN, ARcoef else RETURN, FLOAT(ARcoef)
  110.  
  111. END
  112.  
  113.